This document contains the code used to generate the panels in Figures 3 and 4.

The “big b_a” and “small b_a” parameters identified here were used to generate trajectories for the validation analysis.

The combined figure for Figure 3 was plotted in Matlab using the data generated in “R_0_Surface_Plot_data_for_revision.csv”. The Matlab code is included in the file “Surface_Plot_Panel_Script_add_colorbar_2D.m” in the repo.

library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
source("load_libraries_essential.R")
## Warning: package 'dplyr' was built under R version 3.5.2
library(pomp)
## Warning: package 'pomp' was built under R version 3.5.2
## Welcome to pomp version 2!
## For information on upgrading your pomp version < 2 code, see the
## 'pomp version 2 upgrade guide' at https://kingaa.github.io/pomp/.
source("rahul_theme.R")

Surface plot of b_a vs b_p vs R_0, color by b_a

model_name = "N_12"
load(file = paste0("../Generated_Data/Profiles/", model_name, "_Model/top_2_LL_data/b_a_profile_antibody_surface_plot_data.RData"))

load(file = paste0("../Generated_Data/Profiles/", model_name, "_Model/top_2_LL_data/b_a_profile_antibody_surface_plot_data_with_outlier.RData"))
#head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier)

antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers =
  antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier %>%
  filter(b_p > 0.02)

f1 <- list(
  family = "Arial, sans-serif",
  size = 18,
  color = "black"
)
f2 <- list(
  family = "Old Standard TT, serif",
  size = 15,
  color = "black"
)
R_0_surf_plot_data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%dplyr::select(R_0, b_a, b_p)
write.csv(R_0_surf_plot_data,
          "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/R_0_Surface_Plot_Data_for_Revision.csv",
          row.names = FALSE)

fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
               x = ~b_a, y = ~b_p, z = ~R_0, color = ~R_0)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
                                                zeroline = TRUE,
                                                range = c(0,1),
                                                zerolinecolor = toRGB("black"),
                                                titlefont = f1,
                                                tickfont = f2,
                                  zerolinewidth = 6),
                     yaxis = list(title = ' b_p',
                                  range = c(0,1),
                                  zeroline = TRUE,
                                  showline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6),
                     zaxis = list(title = 'R_0 ',
                                  range = c(0,8.2),
                                  zeroline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6)))

fig
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
                                                zeroline = TRUE,
                                                zerolinecolor = toRGB("black"),
                                                titlefont = f1,
                                                tickfont = f2,
                                  zerolinewidth = 6),
                     yaxis = list(title = ' b_p',
                                  zeroline = TRUE,
                                  showline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6),
                     zaxis = list(title = 'R_0 ',
                                  zeroline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6)))

fig

Surface plot of b_a vs b_p vs R_0_NGM, color by b_a

R_0_NGM_surf_plot_data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%dplyr::select(R_0_NGM, b_a, b_p)
write.csv(R_0_NGM_surf_plot_data,
          "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/R_0_NGM_Surface_Plot_Data_for_Revision.csv",
          row.names = FALSE)

fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
               x = ~b_a, y = ~b_p, z = ~R_0_NGM, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
                                                zeroline = TRUE,
                                                range = c(0,1),
                                                zerolinecolor = toRGB("black"),
                                                titlefont = f1,
                                                tickfont = f2,
                                  zerolinewidth = 6),
                     yaxis = list(title = ' b_p',
                                  range = c(0,1),
                                  zeroline = TRUE,
                                  showline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6),
                     zaxis = list(title = 'R_0_NGM ',
                                  range = c(0,5),
                                  zeroline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6)))

fig
fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
               x = ~b_a, y = ~b_p, z = ~R_0_NGM, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a'),
                     yaxis = list(title = ' b_p'),
                     zaxis = list(title = 'R_0_NGM ')))

fig

Surface plot of R_0 vs R_0_NGM vs b_p, color by b_a

fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
               x = ~R_0, y = ~R_0_NGM, z = ~b_p, color = ~b_a)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' R_0'),
                     yaxis = list(title = ' R_0_NGM'),
                     zaxis = list(title = 'b_p ')))

fig
fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
               x = ~b_a, y = ~b_p, z = ~R_0, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
                                                zeroline = TRUE,
                                                range = c(0,1),
                                                zerolinecolor = toRGB("black"),
                                                titlefont = f1,
                                                tickfont = f2,
                                  zerolinewidth = 6),
                     yaxis = list(title = ' b_p',
                                  range = c(0,1),
                                  zeroline = TRUE,
                                  showline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6),
                     zaxis = list(title = 'R_0_NGM ',
                                  range = c(0,5),
                                  zeroline = TRUE,
                                  zerolinecolor = toRGB("black"),
                                  titlefont = f1,
                                  tickfont = f2,
                                  zerolinewidth = 6)))
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a)) + rahul_man_figure_theme +
  scale_color_viridis_c(name = "Strength \n  of Asymp\n  Transmission") +
  geom_point() +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") 
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Isolate MLE for each scenario

head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)
##          b_a M_0 V_0 K_0      R_0       b_q       b_p       p_S p_H_cond_S
## 1 0.00000000   5  13  14 7.061999 0.2401933 0.7957984 0.1574275  0.1749811
## 2 0.00000000   5  13  14 6.215237 0.2395227 0.9233632 0.1740948  0.1532892
## 3 0.06896552   5  13  14 5.984521 0.2204782 0.8913949 0.1735879  0.1757053
## 4 0.06896552   5  13  14 6.098822 0.2343746 0.9402490 0.1488610  0.1569547
## 5 0.10344828   5  13  14 6.845525 0.2237219 0.7932278 0.1296082  0.1818042
## 6 0.13793103   5  13  14 5.905499 0.2168981 0.8785885 0.1494229  0.1631618
##   phi_E phi_U phi_S   h_V     gamma   N_0      E_0       z_0 C_0
## 1  1.09  1.09   0.2 0.125 58.078038 8e+06 65553.63 10267.580   0
## 2  1.09  1.09   0.2 0.125  8.978508 8e+06 59637.11 11293.462   0
## 3  1.09  1.09   0.2 0.125 35.301204 8e+06 56854.62  9411.456   0
## 4  1.09  1.09   0.2 0.125  6.333218 8e+06 63566.34 13443.034   0
## 5  1.09  1.09   0.2 0.125 48.705248 8e+06 71702.50 12381.941   0
## 6  1.09  1.09   0.2 0.125 17.720241 8e+06 63878.93 12506.505   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2764235
## 2                           17                    22      0.9 0.2750787
## 3                           17                    22      0.9 0.2768054
## 4                           17                    22      0.9 0.2709327
## 5                           17                    22      0.9 0.2752184
## 6                           17                    22      0.9 0.2751406
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        1          10 -627.3202    NA        NA
## 2         0.162 mif1        1          11 -627.2093    NA        NA
## 3         0.162 mif1        2          43 -628.1031    NA        NA
## 4         0.162 mif1        2          44 -627.4300    NA        NA
## 5         0.162 mif1        1          60 -628.1314    NA        NA
## 6         0.162 mif1        2          76 -628.5305    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760        -24.75616    0.004776843            0.1933564
## 2 0.010463251        -26.12591    0.007764530            0.1749998
## 3 0.007808338        -26.08064    0.008118650            0.1775874
## 4 0.008839892        -24.41689    0.002390176            0.1998329
## 5 0.010838884        -25.56156    0.006393602            0.2380878
## 6 0.008314098        -24.34239    0.002107043            0.2090317
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1         3                1                  5         0.01721821
## 2         4                1                  5         0.11137708
## 3        11                1                  5         0.02832765
## 4        13                1                  5         0.15789761
## 5        15                1                  5         0.02053167
## 6        19                1                  5         0.05643264
##   duration_of_symp gamma_total     Beta     R_0_P     R_0_A   R_0_S_1
## 1         5.017218   0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2         5.111377   0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3         5.028328   0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4         5.157898   0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5         5.020532   0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6         5.056433   0.1977679 1.167918 0.9413940 0.6851065 0.8725683
##       R_0_S_2  R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310
range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$E_0)
## [1] 44295.44 71702.50
range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$z_0)
## [1]  9070.932 17695.218
range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_p)
## [1] 0.1400547 0.9865475
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = b_p)) +
  geom_histogram() + rahul_theme
p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

large_b_p_params = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
  filter(b_p >0.90)

p = ggplot(data = large_b_p_params,
         aes(x = R_0,
             y = R_0_NGM,
             color = b_a)) +
  geom_point() + scale_color_viridis_c(name = "Asymp \n Trans \n Strength") +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number \n (Includes Pre-Symp and Asymp Transmission)")  +
  rahul_man_figure_theme +
  theme_white_background+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_high_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p

small_b_p_params = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
  filter(b_p <0.50)
small_b_p_params
##         b_a M_0 V_0 K_0      R_0       b_q       b_p       p_S p_H_cond_S
## 1 0.2413793   5  13  14 8.125304 0.1611435 0.3220870 0.1574198  0.1567379
## 2 0.2758621   5  13  14 7.434055 0.1584788 0.3183492 0.1713446  0.1567255
## 3 0.3448276   5  13  14 6.325641 0.1615416 0.4648878 0.1623848  0.1546026
## 4 0.3448276   5  13  14 6.619633 0.1668452 0.4654342 0.1359111  0.1505722
## 5 0.5517241   5  13  14 5.776447 0.1487695 0.2588775 0.1517047  0.1843260
## 6 0.5862069   5  13  14 6.171486 0.1337064 0.1400547 0.1645830  0.1660020
## 7 0.6896552   5  13  14 4.930404 0.1541091 0.4620370 0.1486280  0.1836424
##   phi_E phi_U phi_S   h_V       gamma   N_0      E_0      z_0 C_0
## 1  1.09  1.09   0.2 0.125   79.534133 8e+06 52522.51 12908.19   0
## 2  1.09  1.09   0.2 0.125 3861.955611 8e+06 50636.94 11633.76   0
## 3  1.09  1.09   0.2 0.125  253.068498 8e+06 46679.50 14443.09   0
## 4  1.09  1.09   0.2 0.125  173.699820 8e+06 53744.73 17695.22   0
## 5  1.09  1.09   0.2 0.125  194.903437 8e+06 56052.18 11154.96   0
## 6  1.09  1.09   0.2 0.125  182.122463 8e+06 44295.44 12335.03   0
## 7  1.09  1.09   0.2 0.125    2.663607 8e+06 55541.02 11957.99   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2784369
## 2                           17                    22      0.9 0.2790559
## 3                           17                    22      0.9 0.2776869
## 4                           17                    22      0.9 0.2808862
## 5                           17                    22      0.9 0.2810877
## 6                           17                    22      0.9 0.2804963
## 7                           17                    22      0.9 0.2827579
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 7 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        2         121 -628.6123    NA        NA
## 2         0.162 mif1        2         138 -628.5431    NA        NA
## 3         0.162 mif1        2         168 -629.0010    NA        NA
## 4         0.162 mif1        1         180 -629.0174    NA        NA
## 5         0.162 mif1        2         271 -628.5576    NA        NA
## 6         0.162 mif1        2         288 -629.1750    NA        NA
## 7         0.162 mif1        2         343 -629.0228    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008142806        -24.85219    0.003761241            0.1912882
## 2 0.009305056        -26.00451    0.006860371            0.1764086
## 3 0.010589918        -25.15982    0.003990094            0.1851563
## 4 0.006043226        -24.79345    0.004001890            0.2221081
## 5 0.004688799        -24.69044    0.002171494            0.1986929
## 6 0.007497558        -25.45496    0.004353770            0.1833871
## 7 0.006280933        -24.50344    0.002076568            0.2082881
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1        31                1                  5       0.0125732181
## 2        36                1                  5       0.0002589362
## 3        43                1                  5       0.0039514993
## 4        46                1                  5       0.0057570583
## 5        68                1                  5       0.0051307458
## 6        70                1                  5       0.0054908109
## 7        76                1                  5       0.3754308107
##   duration_of_symp gamma_total      Beta     R_0_P    R_0_A   R_0_S_1
## 1         5.012573   0.1994983 1.6209846 0.4789891 1.648391 1.2758754
## 2         5.000259   0.1999896 1.4867341 0.4342207 1.699297 1.2737192
## 3         5.003951   0.1998421 1.2641291 0.5391543 1.825610 1.0263769
## 4         5.005757   0.1997700 1.3224039 0.5646715 1.970129 0.8986468
## 5         5.005131   0.1997950 1.1541052 0.2741026 2.700750 0.8754160
## 6         5.005491   0.1997806 1.2329432 0.1584215 3.019029 1.0146071
## 7         5.375431   0.1860316 0.9172109 0.3887939 2.692716 0.6816160
##        R_0_S_2  R_0_NGM
## 1 0.0027054986 3.405961
## 2 0.0000556244 3.407292
## 3 0.0006857403 3.391827
## 4 0.0008789135 3.434326
## 5 0.0007327260 3.851002
## 6 0.0009292432 4.192987
## 7 0.0417811232 3.804907
range(small_b_p_params$b_a)
## [1] 0.2413793 0.6896552
small_b_p_big_b_a = small_b_p_params %>%
  filter(b_a == max(b_a))
small_b_p_small_b_a = small_b_p_params %>%
  filter(b_a == min(b_a))
small_b_p_big_b_a$b_p
## [1] 0.462037
small_b_p_small_b_a$b_p
## [1] 0.322087
hist(small_b_p_params$b_p)

hist(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_p)

big_b_a_values = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
  filter(b_a >= 0.52)
range(big_b_a_values$R_0_NGM)
## [1] 3.153313 4.192987
range(big_b_a_values$R_0)
## [1] 2.927915 6.171486
small_b_a_values =antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
  filter(b_a < 0.52)
range(small_b_a_values$R_0_NGM)
## [1] 2.108494 3.434326
range(small_b_a_values$R_0)
## [1] 3.984487 8.125304
range(antibody_top_2_LL_from_b_a_profile_top_2_LL$p_S)
## [1] 0.1292286 0.1740948
mid_b_p = antibody_top_2_LL_from_b_a_profile_top_2_LL%>%
  filter(b_p <0.55) %>%
  filter(b_p > 0.45) 
mid_b_p_big_b_a = mid_b_p %>%
  filter(b_a == max(b_a))
mid_b_p_small_b_a = mid_b_p %>%
  filter(b_a < max(b_a)) %>%
  filter(b_p == min(b_p))

mid_b_p_big_b_a$R_0_NGM
## [1] 3.804907
mid_b_p_small_b_a$R_0_NGM
## [1] 3.391827
mid_b_p_big_b_a$R_0
## [1] 4.930404
mid_b_p_small_b_a$R_0
## [1] 6.325641
write.csv(mid_b_p_big_b_a, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_mid_b_p_big_b_a_parameter_combination_for_simulation.csv")

write.csv(mid_b_p_small_b_a, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_mid_b_p_small_b_a_parameter_combination_for_simulation.csv")


low_b_p = antibody_top_2_LL_from_b_a_profile_top_2_LL%>%
  filter(b_p < 0.70) %>%
  filter(b_p > 0.50)
hist(low_b_p$b_p)

low_b_p$b_a
## [1] 0.3448276 0.4827586

Pick big b_a parameter combination

big_b_a_MLE = large_b_p_params %>%
  filter(b_a == max(b_a)) %>%
  filter(b_p == max(b_p)) 
big_b_a_MLE
##         b_a M_0 V_0 K_0      R_0      b_q       b_p      p_S p_H_cond_S
## 1 0.9655172   5  13  14 3.083236 0.163873 0.9865475 0.154388  0.1735144
##   phi_E phi_U phi_S   h_V    gamma   N_0      E_0      z_0 C_0
## 1  1.09  1.09   0.2 0.125 11.72791 8e+06 54806.41 11625.01   0
##   social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1                           17                    22      0.9 0.27583
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        1         478 -628.7638    NA        NA
##   loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.01206257        -24.50512    0.001865621            0.2021783
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1        94                1                  5         0.08526666
##   duration_of_symp gamma_total      Beta     R_0_P    R_0_A   R_0_S_1
## 1         5.085267   0.1966465 0.6063076 0.5487626 2.475108 0.4680332
##       R_0_S_2  R_0_NGM
## 1 0.006596615 3.498501
write.csv(big_b_a_MLE, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_big_b_a_parameter_combination_for_simulation.csv")
big_b_a_MLE_comb =   big_b_a_MLE %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)

mid_b_p_big_b_a_comb =   mid_b_p_big_b_a %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)

Load covaraites and obs data

#Load Observed NYC case data
Observed_data = read.csv(paste0(
  "../Generated_Data/observed_data_",
  model_name, ".csv"))
head(Observed_data)
##   Y times obs_prop_positive
## 1 0     1                NA
## 2 0     2        0.00000000
## 3 2     3        0.25000000
## 4 2     4        0.05555556
## 5 7     5        0.15555556
## 6 0     6        0.00000000
### Define start date
true_start_date = as.Date("2020-03-01")
t0 = 0
start_of_year = as.Date("2020-01-01")
first_saturday_in_year = as.Date("2020-01-04")

## Compartment/Queue Cohort Numbers
M = 5
V = 13
K = 14


#Declare Csnippets and data
source("Csnippet_nyc_coronavirus_model_N_12.R")

## Load NYC covariate data
covariate_df = read.csv(file =
                          paste0("../Generated_Data/covariate_data_",
                                 model_name, ".csv"))



### Create covariate table
covar=covariate_table(
  time=covariate_df$times,
  L_advanced_2_days=covariate_df$L_advanced_2_days,
  F_w_y = covariate_df$F_w_y,
  L_orig = covariate_df$L_orig,
  w = covariate_df$Week,
  y = covariate_df$Year,
  times="time"
)

Simulate from big b_a param combination

##Simulation from ML
sim_data_big_b_a = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = big_b_a_MLE_comb,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
#head(sim_data)

sim_data_big_a_FOI_data = sim_data_big_b_a %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
sim_data_big_a_FOI_data = sim_data_big_a_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*big_b_a_MLE$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*big_b_a_MLE$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)

median_FOI_data_big_b_a = sim_data_big_a_FOI_data %>%
  group_by(time) %>%
  summarize(Symptomatic = median(lambda_FOI_S),
            Asymptomatic = median(lambda_FOI_A),
            Presymptomatic = median(lambda_FOI_P)) %>%
  as.data.frame() 
median_FOI_data_big_b_a_melt = median_FOI_data_big_b_a %>%
  melt(id.vars= c("time")) %>%
  dplyr::select(time, FOI = value,
                Infection_Status = variable)

#head(median_FOI_data)
p = ggplot(data = median_FOI_data_big_b_a_melt,
           aes(x = time,
               y = FOI,
               fill = Infection_Status)) +
  geom_area() + rahul_man_figure_theme +
  theme_white_background +
  xlab("Days since March 1, 2020") +
  ylab("Force of Infection") +
  theme(legend.position = c(0.70,0.75)) +
  scale_fill_discrete(name = "Infection Status")+
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Simulate from mid b_p big b_a param combination

##Simulation from ML
sim_data_mid_b_p_big_b_a = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = mid_b_p_big_b_a_comb,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
#head(sim_data)

sim_data_mid_b_p_big_a_FOI_data = sim_data_mid_b_p_big_b_a %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
sim_data_mid_b_p_big_a_FOI_data = sim_data_mid_b_p_big_a_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*mid_b_p_big_b_a_comb$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*mid_b_p_big_b_a_comb$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)

median_FOI_data_mid_b_p_big_b_a = sim_data_mid_b_p_big_a_FOI_data %>%
  group_by(time) %>%
  summarize(Symptomatic = median(lambda_FOI_S),
            Asymptomatic = median(lambda_FOI_A),
            Presymptomatic = median(lambda_FOI_P)) %>%
  as.data.frame() 
median_FOI_data_mid_b_p_big_b_a_melt = median_FOI_data_mid_b_p_big_b_a %>%
  melt(id.vars= c("time")) %>%
  dplyr::select(time, FOI = value,
                Infection_Status = variable)

#head(median_FOI_data)
p = ggplot(data = median_FOI_data_mid_b_p_big_b_a_melt,
           aes(x = time,
               y = FOI,
               fill = Infection_Status)) +
  geom_area() + rahul_man_figure_theme +
  theme_white_background +
  xlab("Days since March 1, 2020") +
  ylab("Force of Infection") +
  theme(legend.position = c(0.70,0.75)) +
  scale_fill_discrete(name = "Infection Status")+
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Identify small b_a param

small_b_a_MLE = large_b_p_params %>%
  filter(b_a > 0) %>%
  filter(b_a == min(b_a)) %>%
  filter(b_p == max(b_p)) 
small_b_a_MLE
##          b_a M_0 V_0 K_0      R_0       b_q      b_p      p_S p_H_cond_S
## 1 0.06896552   5  13  14 6.098822 0.2343746 0.940249 0.148861  0.1569547
##   phi_E phi_U phi_S   h_V    gamma   N_0      E_0      z_0 C_0
## 1  1.09  1.09   0.2 0.125 6.333218 8e+06 63566.34 13443.03   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2709327
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index  loglik nfail trace_num
## 1         0.162 mif1        2          44 -627.43    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008839892        -24.41689    0.002390176            0.1998329
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1        13                1                  5          0.1578976
##   duration_of_symp gamma_total     Beta    R_0_P    R_0_A  R_0_S_1
## 1         5.157898   0.1938774 1.182424 1.019975 0.347037 0.880084
##      R_0_S_2  R_0_NGM
## 1 0.02343045 2.270527
write.csv(small_b_a_MLE, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_small_b_a_parameter_combination_for_simulation.csv")
small_b_a_MLE_comb =   small_b_a_MLE %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)

Identify mid b_p small b_a param

mid_b_p_small_b_a_comb =   mid_b_p_small_b_a %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)

Simulate from small b_a param combination

##Simulation from ML
sim_data_small_b_a = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = small_b_a_MLE_comb,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
#head(sim_data)

sim_data_small_a_FOI_data = sim_data_small_b_a %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
sim_data_small_a_FOI_data = sim_data_small_a_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*small_b_a_MLE$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*small_b_a_MLE$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)

median_FOI_data_small_b_a = sim_data_small_a_FOI_data %>%
  group_by(time) %>%
  summarize(Symptomatic = median(lambda_FOI_S),
            Asymptomatic = median(lambda_FOI_A),
            Presymptomatic = median(lambda_FOI_P)) %>%
  as.data.frame() 
median_FOI_data_small_b_a_melt = median_FOI_data_small_b_a %>%
  melt(id.vars= c("time")) %>%
  dplyr::select(time, FOI = value,
                Infection_Status = variable)

#head(median_FOI_data)
p = ggplot(data = median_FOI_data_small_b_a_melt,
           aes(x = time,
               y = FOI,
               fill = Infection_Status)) +
  geom_area() + rahul_man_figure_theme +
  theme_white_background +
  xlab("Days since March 1, 2020") +
  ylab("Force of Infection") +
  theme(legend.position = c(0.70,0.75)) +
  scale_fill_discrete(name = "Infection Status")+
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Simulate from small b_a param combination

##Simulation from ML
sim_data_mid_b_p_small_b_a = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = mid_b_p_small_b_a_comb,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
#head(sim_data)

sim_data_mid_b_p_small_a_FOI_data = sim_data_mid_b_p_small_b_a %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
sim_data_mid_b_p_small_a_FOI_data = sim_data_mid_b_p_small_a_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*mid_b_p_small_b_a_comb$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*mid_b_p_small_b_a_comb$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)

median_FOI_data_mid_b_p_small_b_a = sim_data_mid_b_p_small_a_FOI_data %>%
  group_by(time) %>%
  summarize(Symptomatic = median(lambda_FOI_S),
            Asymptomatic = median(lambda_FOI_A),
            Presymptomatic = median(lambda_FOI_P)) %>%
  as.data.frame() 
median_FOI_data_mid_b_p_small_b_a_melt = median_FOI_data_mid_b_p_small_b_a %>%
  melt(id.vars= c("time")) %>%
  dplyr::select(time, FOI = value,
                Infection_Status = variable)

#head(median_FOI_data)
p = ggplot(data = median_FOI_data_mid_b_p_small_b_a_melt,
           aes(x = time,
               y = FOI,
               fill = Infection_Status)) +
  geom_area() + rahul_man_figure_theme +
  theme_white_background +
  xlab("Days since March 1, 2020") +
  ylab("Force of Infection") +
  theme(legend.position = c(0.70,0.75)) +
  scale_fill_discrete(name = "Infection Status")+
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.jpeg",
     quality = 100)
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
outlier = antibody_top_2_LL_from_b_a_profile_top_2_LL %>%
  filter(R_0 > 15)
outlier
##          b_a M_0 V_0 K_0      R_0      b_q         b_p       p_S
## 1 0.06896552   5  13  14 17.77109 0.133977 0.008385475 0.1630488
##   p_H_cond_S phi_E phi_U phi_S   h_V    gamma   N_0      E_0      z_0 C_0
## 1  0.1729462  1.09  1.09   0.2 0.125 960.6433 8e+06 54219.45 11091.14   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2807117
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        1          37 -629.0826    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.004121517        -25.32007    0.007746291            0.1847238
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1        10                1                  5        0.001040969
##   duration_of_symp gamma_total     Beta      R_0_P    R_0_A  R_0_S_1
## 1         5.001041   0.1999584 3.553478 0.02733725 1.025547 2.896952
##        R_0_S_2  R_0_NGM
## 1 0.0004988189 3.950335

Plot b_a * b_p

product_data =
  antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
  mutate(b_a_times_b_p = b_a*b_p,
         b_a_plus_b_p = b_a + b_p)
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a_times_b_p)) + rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) +
   theme(axis.title.x = element_text(face = "plain",
                                    size = 24),
        axis.title.y = element_text(face = "plain",
                                    size = 24))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a_times_b_p)) + rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) +
   theme(axis.title.x = element_text(face = "plain",
                                    size = 24),
        axis.title.y = element_text(face = "plain",
                                    size = 24))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a_plus_b_p)) + rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point() +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) +
   theme(axis.title.x = element_text(face = "plain",
                                    size = 24),
        axis.title.y = element_text(face = "plain",
                                    size = 24))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_plus_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a_plus_b_p)) + rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point() +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_plus_b_p_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a)) +
  rahul_man_figure_theme +
  scale_color_viridis_c(name = expression(b[a])) +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  theme_white_background +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number")  +
  theme(axis.title.x = element_text(face = "plain",
                                    size = 24),
        axis.title.y = element_text(face = "plain",
                                    size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))#+ #+
  # theme(legend.position = c(0.75,0.75)) 
 # labs(y = expression(R["0 NGM"]),
 #     x = expression(R[0]))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a)) +
  rahul_man_figure_theme +
  scale_color_viridis_c(name = expression(b[a])) +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  theme_white_background +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number")  +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))+ #+
  # theme(legend.position = c(0.75,0.75)) 
  labs(y = expression(R[0["NGM"]]),
       x = expression(R[0]))+
  theme(axis.title.x = element_text(face = "plain",
                                    size = 32),
        axis.title.y = element_text(face = "plain",
                                    size = 32))
p

#expression(R[0["NGM"]])

png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_symb_plot.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
write.csv(product_data,
          "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/product_data_for_Fig_3C.csv", row.names = FALSE)

p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a)) +
  rahul_man_figure_theme +
  scale_color_viridis_c(name = expression(b[a])) +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  theme_white_background +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number")  +
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) #+
  # theme(legend.position = c(0.75,0.75))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_p)) +
  rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point() +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_p)) +
  rahul_man_figure_theme +
  scale_color_viridis_c() +
  geom_point() +
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_p_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a*b_p)) +
  rahul_man_figure_theme +
  scale_color_viridis_c(name = expression(paste(b[a], "*", b[p]))) +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24)) +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_a*b_p)) +
  rahul_man_figure_theme +
  scale_color_viridis_c(name = expression(paste(b[a], "*", b[p]))) +
  geom_point(size = 5) +
  scale_x_continuous(breaks=c(seq(2,9))) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, 
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) + 
  xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
  theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 24),
        axis.title.y = element_text(face = "plain", size = 24))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) #+
  # theme(legend.position = c(0.75, 0.75))
p

png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = b_a*b_p,
               y = R_0)) +
  rahul_man_figure_theme +
  geom_point(size= 5) +
  ylab("Secondary cases from each \n primary symptomatic case") +
  xlab(expression(paste(b[a],'*',b[p]))) +
  theme_white_background+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_times_b_p_vs_R_0_antibody_LL_top_no_outliers.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = b_a*b_p,
               y = R_0)) +
  rahul_man_figure_theme +
  geom_point(size= 5) +
  ylab("Secondary cases from each \n primary symptomatic case") +
  xlab(expression(paste(b[a],'*',b[p]))) +
  theme_white_background+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_times_b_p_vs_R_0_antibody_LL_top_no_outliers_circle_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = b_a+b_p,
               y = R_0)) +
  rahul_man_figure_theme +
  geom_point() +
  ylab("Secondary cases from each \n primary symptomatic case") +
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_plus_b_p_vs_R_0_antibody_LL_top_no_outliers.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = product_data,
           aes(x = b_a+b_p,
               y = R_0)) +
  rahul_man_figure_theme +
  geom_point() +
  ylab("Secondary cases from each \n primary symptomatic case")+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
  theme(axis.text.x = element_text(size=21)) +
  theme(axis.text.y = element_text(size=21)) 
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_plus_b_p_vs_R_0_antibody_LL_top_no_outliers_cirlce_mid_b_p.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)
##          b_a M_0 V_0 K_0      R_0       b_q       b_p       p_S p_H_cond_S
## 1 0.00000000   5  13  14 7.061999 0.2401933 0.7957984 0.1574275  0.1749811
## 2 0.00000000   5  13  14 6.215237 0.2395227 0.9233632 0.1740948  0.1532892
## 3 0.06896552   5  13  14 5.984521 0.2204782 0.8913949 0.1735879  0.1757053
## 4 0.06896552   5  13  14 6.098822 0.2343746 0.9402490 0.1488610  0.1569547
## 5 0.10344828   5  13  14 6.845525 0.2237219 0.7932278 0.1296082  0.1818042
## 6 0.13793103   5  13  14 5.905499 0.2168981 0.8785885 0.1494229  0.1631618
##   phi_E phi_U phi_S   h_V     gamma   N_0      E_0       z_0 C_0
## 1  1.09  1.09   0.2 0.125 58.078038 8e+06 65553.63 10267.580   0
## 2  1.09  1.09   0.2 0.125  8.978508 8e+06 59637.11 11293.462   0
## 3  1.09  1.09   0.2 0.125 35.301204 8e+06 56854.62  9411.456   0
## 4  1.09  1.09   0.2 0.125  6.333218 8e+06 63566.34 13443.034   0
## 5  1.09  1.09   0.2 0.125 48.705248 8e+06 71702.50 12381.941   0
## 6  1.09  1.09   0.2 0.125 17.720241 8e+06 63878.93 12506.505   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2764235
## 2                           17                    22      0.9 0.2750787
## 3                           17                    22      0.9 0.2768054
## 4                           17                    22      0.9 0.2709327
## 5                           17                    22      0.9 0.2752184
## 6                           17                    22      0.9 0.2751406
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        1          10 -627.3202    NA        NA
## 2         0.162 mif1        1          11 -627.2093    NA        NA
## 3         0.162 mif1        2          43 -628.1031    NA        NA
## 4         0.162 mif1        2          44 -627.4300    NA        NA
## 5         0.162 mif1        1          60 -628.1314    NA        NA
## 6         0.162 mif1        2          76 -628.5305    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760        -24.75616    0.004776843            0.1933564
## 2 0.010463251        -26.12591    0.007764530            0.1749998
## 3 0.007808338        -26.08064    0.008118650            0.1775874
## 4 0.008839892        -24.41689    0.002390176            0.1998329
## 5 0.010838884        -25.56156    0.006393602            0.2380878
## 6 0.008314098        -24.34239    0.002107043            0.2090317
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1         3                1                  5         0.01721821
## 2         4                1                  5         0.11137708
## 3        11                1                  5         0.02832765
## 4        13                1                  5         0.15789761
## 5        15                1                  5         0.02053167
## 6        19                1                  5         0.05643264
##   duration_of_symp gamma_total     Beta     R_0_P     R_0_A   R_0_S_1
## 1         5.017218   0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2         5.111377   0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3         5.028328   0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4         5.157898   0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5         5.020532   0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6         5.056433   0.1977679 1.167918 0.9413940 0.6851065 0.8725683
##       R_0_S_2  R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = p_H_cond_S,  y = b_a)) +
  geom_point() + rahul_man_figure_theme
p

p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = p_S,  y = b_a)) +
  geom_point() + rahul_man_figure_theme
p

head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)
##          b_a M_0 V_0 K_0      R_0       b_q       b_p       p_S p_H_cond_S
## 1 0.00000000   5  13  14 7.061999 0.2401933 0.7957984 0.1574275  0.1749811
## 2 0.00000000   5  13  14 6.215237 0.2395227 0.9233632 0.1740948  0.1532892
## 3 0.06896552   5  13  14 5.984521 0.2204782 0.8913949 0.1735879  0.1757053
## 4 0.06896552   5  13  14 6.098822 0.2343746 0.9402490 0.1488610  0.1569547
## 5 0.10344828   5  13  14 6.845525 0.2237219 0.7932278 0.1296082  0.1818042
## 6 0.13793103   5  13  14 5.905499 0.2168981 0.8785885 0.1494229  0.1631618
##   phi_E phi_U phi_S   h_V     gamma   N_0      E_0       z_0 C_0
## 1  1.09  1.09   0.2 0.125 58.078038 8e+06 65553.63 10267.580   0
## 2  1.09  1.09   0.2 0.125  8.978508 8e+06 59637.11 11293.462   0
## 3  1.09  1.09   0.2 0.125 35.301204 8e+06 56854.62  9411.456   0
## 4  1.09  1.09   0.2 0.125  6.333218 8e+06 63566.34 13443.034   0
## 5  1.09  1.09   0.2 0.125 48.705248 8e+06 71702.50 12381.941   0
## 6  1.09  1.09   0.2 0.125 17.720241 8e+06 63878.93 12506.505   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2764235
## 2                           17                    22      0.9 0.2750787
## 3                           17                    22      0.9 0.2768054
## 4                           17                    22      0.9 0.2709327
## 5                           17                    22      0.9 0.2752184
## 6                           17                    22      0.9 0.2751406
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index    loglik nfail trace_num
## 1         0.162 mif1        1          10 -627.3202    NA        NA
## 2         0.162 mif1        1          11 -627.2093    NA        NA
## 3         0.162 mif1        2          43 -628.1031    NA        NA
## 4         0.162 mif1        2          44 -627.4300    NA        NA
## 5         0.162 mif1        1          60 -628.1314    NA        NA
## 6         0.162 mif1        2          76 -628.5305    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760        -24.75616    0.004776843            0.1933564
## 2 0.010463251        -26.12591    0.007764530            0.1749998
## 3 0.007808338        -26.08064    0.008118650            0.1775874
## 4 0.008839892        -24.41689    0.002390176            0.1998329
## 5 0.010838884        -25.56156    0.006393602            0.2380878
## 6 0.008314098        -24.34239    0.002107043            0.2090317
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1         3                1                  5         0.01721821
## 2         4                1                  5         0.11137708
## 3        11                1                  5         0.02832765
## 4        13                1                  5         0.15789761
## 5        15                1                  5         0.02053167
## 6        19                1                  5         0.05643264
##   duration_of_symp gamma_total     Beta     R_0_P     R_0_A   R_0_S_1
## 1         5.017218   0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2         5.111377   0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3         5.028328   0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4         5.157898   0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5         5.020532   0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6         5.056433   0.1977679 1.167918 0.9413940 0.6851065 0.8725683
##       R_0_S_2  R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = p_S)) +
  geom_histogram() +
  rahul_man_figure_theme
p
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$p_S)
## [1] 0.1292286 0.1740948

Averaging over time

FOI_df = data.frame(matrix(nrow = 0, ncol =
                             ncol(                              antibody_top_2_LL_from_b_a_profile_top_2_LL) + 7))
colnames(FOI_df) = c(colnames(
  antibody_top_2_LL_from_b_a_profile_top_2_LL
), "Mean_Symp_prop", "Mean_Asymp_prop", "Mean_Presymp_Prop", "Mean_Symp", "Mean_Asymp", "Mean_Presymp", "Average_Total_FOI")
for(param_index in seq(1:nrow(antibody_top_2_LL_from_b_a_profile_top_2_LL))){
  print(param_index)
  param_combo =
    antibody_top_2_LL_from_b_a_profile_top_2_LL[param_index,]
  param_combo_for_sim = param_combo %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)
  
  sim_data_big_single_param = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = param_combo_for_sim,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
#head(sim_data)

sim_data_single_param_FOI_data = sim_data_big_single_param %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)

median_FOI_data_single_param = sim_data_single_param_FOI_data %>%
  group_by(time) %>%
  summarize(Symptomatic = median(lambda_FOI_S),
            Asymptomatic = median(lambda_FOI_A),
            Presymptomatic = median(lambda_FOI_P)) %>%
  as.data.frame() %>%
  mutate(Total = Symptomatic + Asymptomatic + Presymptomatic) %>%
  mutate(Symp_Prop = Symptomatic/Total,
         Asymp_Prop = Asymptomatic/Total,
         Presymp_Prop = Presymptomatic/Total)

single_param_avg_FOI = median_FOI_data_single_param %>%
  summarize(Mean_Symp_prop = mean(Symp_Prop),
            Mean_Asymp_Prop = mean(Asymp_Prop),
            Mean_Presymp_Prop = mean(Presymp_Prop),
            Mean_Symp = mean(Symptomatic),
            Mean_Asymp = mean(Asymptomatic),
            Mean_Presymp = mean(Presymptomatic),
            Average_Total_FOI = mean(Total))
single_param_FOI_df = cbind(param_combo, single_param_avg_FOI)
FOI_df = rbind(FOI_df,single_param_FOI_df)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36
p = ggplot(data = FOI_df,
           aes(x = b_a, y = b_p,
               color = Mean_Asymp_Prop)) +
  geom_point() + scale_color_viridis_c() +
  rahul_man_figure_theme + theme_white_background
p

p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Symp)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Symp_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Presymp)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Presymp_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Asymp)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Asymp_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = Average_Total_FOI)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Average_Total_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = b_p,
               color = Mean_Symp_prop)) +
  geom_point() + scale_color_viridis_c() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_b_p_color_Mean_symp_prop_of_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = b_p,
               color = (R_0_S_1 + R_0_S_2)/R_0_NGM)) +
  geom_point() + scale_color_viridis_c(name = "Prop Symp R_0") +
  rahul_man_figure_theme + theme_white_background
p

p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Symp_prop)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_symp_prop_of_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Presymp_Prop)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Presymp_prop_of_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = Mean_Asymp_Prop)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Asymp_prop_of_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = (R_0_S_1 + R_0_S_2)/R_0_NGM)) +
  geom_point() + ylab("Prop Symp R_0") +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Symp_prop_of_R_0.png")
print(p)
dev.off
## function (which = dev.cur()) 
## {
##     if (which == 1) 
##         stop("cannot shut down device 1 (the null device)")
##     .External(C_devoff, as.integer(which))
##     dev.cur()
## }
## <bytecode: 0x7f955d574698>
## <environment: namespace:grDevices>
p = ggplot(data = FOI_df,
           aes(x = b_a, y = (R_0_P)/R_0_NGM)) +
  geom_point() + ylab("Prop Pre-Symp R_0") +
  rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Pre_Symp_prop_of_R_0.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a, y = (R_0_A)/R_0_NGM)) +
  geom_point() + ylab("Prop Asymp R_0") +
  rahul_man_figure_theme + theme_white_background
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Asymp_prop_of_R_0.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data = FOI_df,
           aes(x = b_a*b_p,
               y = Mean_Symp_prop)) +
  geom_point() +
  rahul_man_figure_theme + theme_white_background
p

Figure 4 Alternate

Sub-funcitons

get_median_and_quantiles = function(single_b_a_value_FOI_data, type){
    if(type == "Symptomatic"){
      single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
        mutate(FOI_contrib = lambda_FOI_S/lambda_FOI_total)
    }else{
        if(type == "Asymptomatic"){
          single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
        mutate(FOI_contrib = lambda_FOI_A/lambda_FOI_total)
        }else{
          if(type == "Presymptomatic"){
            single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
        mutate(FOI_contrib = lambda_FOI_P/lambda_FOI_total)
          }else{
            error("Invalid infection class")
          }
        }
      }
      single_b_a_value_summary = single_b_a_value_FOI_data %>%
      group_by(time) %>%
      summarize(median_FOI_contrib = 
                  median(FOI_contrib),
                upper_quantile_FOI_contrib = quantile(FOI_contrib, probs = 0.975),
                lower_quantile_FOI_contrib = quantile(FOI_contrib, probs = 0.025)
                ) %>%
        as.data.frame()%>%
        mutate(Infection_Type = type)
}

Define sampling times for FOI

head(Observed_data)
##   Y times obs_prop_positive
## 1 0     1                NA
## 2 0     2        0.00000000
## 3 2     3        0.25000000
## 4 2     4        0.05555556
## 5 7     5        0.15555556
## 6 0     6        0.00000000
peak_time = Observed_data %>%
  filter(Y == max(Y)) %>%
  dplyr::select(times) %>%
  as.numeric()

pre_peak_time = peak_time - 28
post_peak_time = peak_time + 28
p = ggplot(data = Observed_data, aes(x = times, y = Y)) + geom_point() + geom_line() + rahul_man_figure_theme +
  geom_vline(xintercept = pre_peak_time, color = 'blue',
             size = 1.5) +
  geom_vline(xintercept = peak_time, color = 'red',
             size = 1.5) +
  geom_vline(xintercept = post_peak_time, color = 'orange',
             size = 1.5) +
  theme_white_background
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/sampling_points_for_FOI.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Calculate summary statistics for FOI at each sampling time for each \(b_a\) including outlier value

b_a_points = unique(antibody_top_2_LL_from_b_a_profile_top_2_LL$b_a)

FOI_summary = data.frame(matrix(nrow = 0, ncol =6 ))
colnames(FOI_summary) = c("time" ,"median_FOI_contrib",
                          "upper_quantile_FOI_contrib",
                          "lower_quantile_FOI_contrib",
                          "Infection_Type", "b_a") 
for(b_a_index in seq(1:length(b_a_points))){
  print(b_a_index)
  single_b_a_value = b_a_points[b_a_index]
  param_combs_with_b_a_index_value =
    antibody_top_2_LL_from_b_a_profile_top_2_LL %>%
    filter(b_a == single_b_a_value)
  
  pre_peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(pre_peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  post_peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(post_peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  
  for(param_index in
      seq(1:nrow(param_combs_with_b_a_index_value))){
    single_param_combo =
      param_combs_with_b_a_index_value[param_index,]
    param_combo_for_sim = single_param_combo %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)
  
  sim_data_single_param = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = param_combo_for_sim,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
  #head(sim_data)
  
  sim_data_single_param_FOI_data = sim_data_single_param %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
  sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P) %>%
  dplyr::select(time, .id, lambda_FOI_S, lambda_FOI_A, lambda_FOI_P,
                lambda_FOI_total) %>%
  mutate(param_index = param_index)
  pre_peak_single_param_FOI_data =
  sim_data_single_param_FOI_data %>%
  filter(time == pre_peak_time)
  peak_time_single_param_FOI_data = 
    sim_data_single_param_FOI_data %>%
    filter(time == peak_time)
  post_peak_single_param_FOI_data =
    sim_data_single_param_FOI_data %>%
    filter(time == post_peak_time)

  pre_peak_single_b_a_value_FOI_data = rbind(
    pre_peak_single_b_a_value_FOI_data,
    pre_peak_single_param_FOI_data
)

peak_single_b_a_value_FOI_data = rbind(
  peak_single_b_a_value_FOI_data,
  peak_time_single_param_FOI_data
)

post_peak_single_b_a_value_FOI_data = rbind(
  post_peak_single_b_a_value_FOI_data,
  post_peak_single_param_FOI_data
)

    
    
  }
  
  ### Calculate % contribution of FOI from each category 
  single_b_a_value_FOI_data = rbind(
    pre_peak_single_b_a_value_FOI_data,
    peak_single_b_a_value_FOI_data) %>%
    rbind(post_peak_single_b_a_value_FOI_data)
  single_b_a_Presymp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Presymptomatic")
  
  single_b_a_Symp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Symptomatic")
  
  single_b_a_Asymp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Asymptomatic")
  
  single_b_a_FOI_summary = rbind(single_b_a_Presymp_summary,
                                 single_b_a_Symp_summary) %>%
    rbind(single_b_a_Asymp_summary) %>%
    mutate(b_a = single_b_a_value)

  FOI_summary = FOI_summary %>%
    rbind(single_b_a_FOI_summary)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21

Isolate times

peak_FOI_summary = FOI_summary %>%
  filter(time == peak_time)

pre_peak_FOI_summary = FOI_summary %>%
  filter(time == pre_peak_time)
post_peak_FOI_summary = FOI_summary %>%
  filter(time == post_peak_time)

Plot FOI

Peak

p = ggplot(data= peak_FOI_summary,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection at Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_peak_include_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Pre-Peak

p = ggplot(data= pre_peak_FOI_summary,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection 4 Weeks Before Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_pre_peak_include_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Post-Peak

p = ggplot(data= post_peak_FOI_summary,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection 4 Weeks After Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_post_peak_include_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Calculate summary statistics for FOI at each sampling time for each \(b_a\) excluding outlier value

second_outlier = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier %>%
  filter(b_p < 0.02)



b_a_points_no_outlier = unique(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_a)

FOI_summary_no_outlier = data.frame(matrix(nrow = 0, ncol =6 ))
colnames(FOI_summary_no_outlier) = c("time" ,"median_FOI_contrib",
                          "upper_quantile_FOI_contrib",
                          "lower_quantile_FOI_contrib",
                          "Infection_Type", "b_a") 
for(b_a_index in seq(1:length(b_a_points_no_outlier))){
  print(b_a_index)
  single_b_a_value = b_a_points_no_outlier[b_a_index]
  param_combs_with_b_a_index_value_no_outlier =
    antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
    filter(b_a == single_b_a_value)
  
  pre_peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(pre_peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  post_peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(post_peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  peak_single_b_a_value_FOI_data =
    data.frame(matrix(nrow = 0, ncol = 13))
  colnames(peak_single_b_a_value_FOI_data) =
    c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
      "lambda_FOI_P","lambda_FOI_total", "param_index")
  
  for(param_index in
      seq(1:nrow(param_combs_with_b_a_index_value_no_outlier))){
    single_param_combo =
      param_combs_with_b_a_index_value_no_outlier[param_index,]
    param_combo_for_sim = single_param_combo %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)
  
  sim_data_single_param = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = param_combo_for_sim,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")
  #head(sim_data)
  
  sim_data_single_param_FOI_data = sim_data_single_param %>%
  dplyr::select(time, .id, I_P, I_S_1, I_S_2,
                I_A, N, beta_t) 
  sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
  mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
         lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
         lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
  mutate(lambda_FOI_total =
           lambda_FOI_S + lambda_FOI_A + lambda_FOI_P) %>%
  dplyr::select(time, .id, lambda_FOI_S, lambda_FOI_A, lambda_FOI_P,
                lambda_FOI_total) %>%
  mutate(param_index = param_index)
  pre_peak_single_param_FOI_data =
  sim_data_single_param_FOI_data %>%
  filter(time == pre_peak_time)
  peak_time_single_param_FOI_data = 
    sim_data_single_param_FOI_data %>%
    filter(time == peak_time)
  post_peak_single_param_FOI_data =
    sim_data_single_param_FOI_data %>%
    filter(time == post_peak_time)

  pre_peak_single_b_a_value_FOI_data = rbind(
    pre_peak_single_b_a_value_FOI_data,
    pre_peak_single_param_FOI_data
)

peak_single_b_a_value_FOI_data = rbind(
  peak_single_b_a_value_FOI_data,
  peak_time_single_param_FOI_data
)

post_peak_single_b_a_value_FOI_data = rbind(
  post_peak_single_b_a_value_FOI_data,
  post_peak_single_param_FOI_data
)

    
    
  }
  
  ### Calculate % contribution of FOI from each category 
  single_b_a_value_FOI_data = rbind(
    pre_peak_single_b_a_value_FOI_data,
    peak_single_b_a_value_FOI_data) %>%
    rbind(post_peak_single_b_a_value_FOI_data)
  single_b_a_Presymp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Presymptomatic")
  
  single_b_a_Symp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Symptomatic")
  
  single_b_a_Asymp_summary = get_median_and_quantiles(
    single_b_a_value_FOI_data = single_b_a_value_FOI_data, 
    type = "Asymptomatic")
  
  single_b_a_FOI_summary = rbind(single_b_a_Presymp_summary,
                                 single_b_a_Symp_summary) %>%
    rbind(single_b_a_Asymp_summary) %>%
    mutate(b_a = single_b_a_value)

  FOI_summary_no_outlier = FOI_summary_no_outlier %>%
    rbind(single_b_a_FOI_summary)
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21

Isolate times

peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
  filter(time == peak_time)

pre_peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
  filter(time == pre_peak_time)
post_peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
  filter(time == post_peak_time)

Peak

p = ggplot(data= peak_FOI_summary_no_outlier,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection at Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_peak_no_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Pre-Peak

p = ggplot(data= pre_peak_FOI_summary_no_outlier,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection 4 Weeks Before Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_pre_peak_no_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Post-Peak

p = ggplot(data= post_peak_FOI_summary_no_outlier,
           aes(x = b_a,
               y = median_FOI_contrib,
               color = Infection_Type)) +
  geom_point() +
  geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
                                ymax = upper_quantile_FOI_contrib,
                                fill = Infection_Type), alpha = 0.5) +
  rahul_man_figure_theme +
  theme_white_background + xlab("Asymptomatic  Transmission Rate ($b_a$) ") +
  ylab("% Contribution to the Force of Infection 4 Weeks After Peak")
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_post_peak_no_outlier.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Peak median density plot

p = ggplot(data=peak_FOI_summary_no_outlier,
       aes(x=factor(b_a), y=median_FOI_contrib, fill=Infection_Type)) +
    geom_bar(stat="identity", position="stack") +
  geom_errorbar(aes(ymin = lower_quantile_FOI_contrib,
                ymax = upper_quantile_FOI_contrib))
p

### Sub-function for stacking data

library(tidyr)
## Warning: package 'tidyr' was built under R version 3.5.2
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
stack_FOI_error_bars = function(FOI_summary_no_outlier){
  
test_df_asymp =FOI_summary_no_outlier %>%
  filter(Infection_Type == "Asymptomatic")

spreaded_median_FOI = FOI_summary_no_outlier %>%
  spread(Infection_Type, median_FOI_contrib) 
median_asymp_FOI = spreaded_median_FOI %>%
  dplyr::select(b_a, Asymptomatic_Median = Asymptomatic) %>%
  na.omit()
median_symp_FOI = spreaded_median_FOI %>%
  dplyr::select(b_a, Symptomatic_Median = Symptomatic) %>%
  na.omit()

median_presymp_FOI = spreaded_median_FOI %>%
  dplyr::select(b_a, Presymptomatic_Median = Presymptomatic) %>%
  na.omit()

median_combined_FOI = median_asymp_FOI %>%
  join(median_presymp_FOI) %>%
  join(median_symp_FOI)

FOI_summary_no_outlier_with_eb = FOI_summary_no_outlier %>%
  join(median_combined_FOI) %>%
  mutate(Add_Presymp = as.numeric(Infection_Type == "Asymptomatic"),
         Add_Symp = as.numeric(Infection_Type %in% c("Asymptomatic", "Presymptomatic" )),
         Stacked_Median_Adjustment = Add_Presymp*Presymptomatic_Median +
           Add_Symp*Symptomatic_Median)

FOI_summary_no_outlier_with_eb = FOI_summary_no_outlier_with_eb%>%
  mutate(b_a_factor = as.factor(b_a))

return(FOI_summary_no_outlier_with_eb)
}

Stack error bars

peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
  FOI_summary_no_outlier = peak_FOI_summary_no_outlier)
## Joining by: b_a
## Joining by: b_a
## Joining by: b_a
pre_peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
  FOI_summary_no_outlier = pre_peak_FOI_summary_no_outlier)
## Joining by: b_a
## Joining by: b_a
## Joining by: b_a
post_peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
  FOI_summary_no_outlier = post_peak_FOI_summary_no_outlier
)
## Joining by: b_a
## Joining by: b_a
## Joining by: b_a
all_b_a_levels = levels(peak_FOI_summary_no_outlier_with_eb$b_a_factor)
b_a_level_tickmarks = all_b_a_levels[c(1,7,13,18,21)]
b_a_tick_labels = c("0", "0.24", "0.52", "0.79", "0.97")

Plot stacked bar plot for peak

p = ggplot(data=peak_FOI_summary_no_outlier_with_eb,
       aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
    geom_bar(stat="identity", position="stack") +
  geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
                ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
  scale_x_discrete(breaks = b_a_level_tickmarks,
                   labels = b_a_tick_labels) +
               xlab('Relative Asymptomatic  \n Transmission Rate') +
  ylab("Contribution to  Force of Infection At Peak ") +
  scale_fill_discrete(name = "Infection \n Type") +
  rahul_man_figure_theme + theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 12),
        axis.title.y = element_text(face = "plain", size = 12),
        axis.text.x = element_text(face = "plain", size = 11),
        axis.text.y = element_text(face = "plain", size = 11),
        legend.text = element_text(face = "plain", size = 11),
        legend.title = element_text(face = "plain", size = 12))+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) 
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Peak_exclude_two_outliers.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2
pdf("../Figures/Profiles/N_12_Model/Man_Figs/stacked_bar_plot_b_a_vs_cont_to_FOI_at_peak_exclude_two_outliers.pdf")
print(p)
dev.off()
## quartz_off_screen 
##                 2
p = ggplot(data=peak_FOI_summary_no_outlier_with_eb,
       aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
    geom_bar(stat="identity", position="stack") +
  geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
                ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
  scale_x_discrete(breaks = b_a_level_tickmarks,
                   labels = b_a_tick_labels) +
               xlab('Relative Asymptomatic  \n Transmission Rate') +
  ylab("Contribution to  Force of Infection At Peak ") +
  scale_fill_discrete(name = "Infection \n Type") +
  rahul_man_figure_theme + theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 12),
        axis.title.y = element_text(face = "plain", size = 12),
        axis.text.x = element_text(face = "plain", size = 11),
        axis.text.y = element_text(face = "plain", size = 11),
        legend.text = element_text(face = "plain", size = 11),
        legend.title = element_text(face = "plain", size = 12)) +
  scale_fill_grey(name = "Infection \n Type")+
  theme(axis.line = element_line(colour = 'black', size = 1))+
  theme(axis.ticks = element_line(colour = "black", size = 1.5)) 
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
p

pdf("../Figures/Profiles/N_12_Model/Man_Figs/stacked_bar_plot_b_a_vs_cont_to_FOI_at_peak_exclude_two_outliers_greyscale.pdf")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Plot stacked bar plot for 4 weeks before peak

p = ggplot(data=pre_peak_FOI_summary_no_outlier_with_eb,
       aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
    geom_bar(stat="identity", position="stack") +
  geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
                ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
  scale_x_discrete(breaks = b_a_level_tickmarks,
                   labels = b_a_tick_labels) +
               xlab('Relative Asymptomatic  \n Transmission Rate') +
  ylab("Contribution to  Force of Infection \n 4 Weeks Before Peak ") +
  scale_fill_discrete(name = "Infection \n Type") +
  rahul_man_figure_theme + theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 12),
        axis.title.y = element_text(face = "plain", size = 12),
        axis.text.x = element_text(face = "plain", size = 11),
        axis.text.y = element_text(face = "plain", size = 11),
        legend.text = element_text(face = "plain", size = 11),
        legend.title = element_text(face = "plain", size = 12))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Pre_Peak_exclude_two_outliers.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Plot stacked bar plot for 4 weeks after peak

p = ggplot(data=post_peak_FOI_summary_no_outlier_with_eb,
       aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
    geom_bar(stat="identity", position="stack") +
  geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
                ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
  scale_x_discrete(breaks = b_a_level_tickmarks,
                   labels = b_a_tick_labels) +
               xlab('Relative Asymptomatic  \n Transmission Rate') +
  ylab("Contribution to  Force of Infection \n 4 Weeks After Peak ") +
  scale_fill_discrete(name = "Infection \n Type") +
  rahul_man_figure_theme + theme_white_background +
  theme(axis.title.x = element_text(face = "plain", size = 12),
        axis.title.y = element_text(face = "plain", size = 12),
        axis.text.x = element_text(face = "plain", size = 11),
        axis.text.y = element_text(face = "plain", size = 11),
        legend.text = element_text(face = "plain", size = 11),
        legend.title = element_text(face = "plain", size = 12))
p

png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Post_Peak_exclude_two_outliers.png")
print(p)
dev.off()
## quartz_off_screen 
##                 2

Sup Figures

p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
           aes(x = R_0,
               y = R_0_NGM,
               color = b_p)) +
  geom_point(size = 5) +
  scale_color_viridis_c() +
  rahul_man_figure_theme +
  theme_white_background +
  scale_x_continuous(breaks=c(seq(2,10,1), 15, 18)) +
  scale_y_continuous(breaks=seq(2,5,1)) +
  coord_cartesian(expand = FALSE, #turn off axis expansion (padding)
                  xlim = c(1.75, 9), ylim = c(1.75, 5.25)) #manually set limits
p

png(paste0("../Figures/Profiles/", model_name, "_Model/Sup_Figs/",
           "R_0_vs_R_0_NGM_color_by_b_p", model_name,
           "_model_top_2_antibody_LL_from_b_a_profile_peak_LL_no_2_outliers.png"))
print(p)
dev.off()
## quartz_off_screen 
##                 2